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Abstract 

Many challenging tasks in sensor networks, including sensor calibration, rank- 
ing of nodes, monitoring, event region detection, collaborative filtering, col- 
laborative signal processing, etc., can be formulated as a problem of solving a 
linear system of equations. Several recent works propose different distributed 
algorithms for solving these problems, usually by using linear iterative nu- 
merical methods. 

The main problem with previous approaches is that once the problem in- 
puts change during the process of computation, the computation may output 
unexpected results. In real life settings, sensor measurements are subject to 
varying environmental conditions and to measurement noise. 

We present a simple iterative scheme called SS- Iterative for solving 
systems of linear equations, and examine its properties in the self- stabilizing 
perspective. We analyze the behavior of the proposed scheme under changing 
input sequences using two different assumptions on the input: a box bound, 
and a probabilistic distribution. 

As study, we discuss the sensor calibration problem and provide 

simulation results to support the applicability of our approach. 
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1. Introduction 

Many challenging tasks in sensor networks, for example distributed rank- 
ing algorithms of nodes and data items [4], collaborative filtering [1], localiza- 
tion [14], collaborative signal processing [16], region detection [13], etc., can 
be formulated as a problem of solving a linear system of equations. Several 
recent works [14], [16], [13] propose different distributed algorithms for solving 
these problems, usually by linear iterative numerical methods. 

In this work, we extend the settings of the above approaches by adding 
another dimension to the problem. Specifically, we are interested in self- 
stabilizing algorithms, that continuously run and converge to a solution from 
any initial state. This aspect of the problem is highly important due to the 
dynamic nature of the network and the frequent changes in the values being 
tracked. 

As study, we show that the calibration of local sensors' readings can 

be formulated as a linear system of equations Ax = b, where x represents the 
calibrated output reading, b represents the local reading, and A represents 
a weighted communication graph. However, our work is general and can be 
applied to any problem that can be formulated as a distributed solution to 
a linear system of equations, including the previous works mentioned above. 

The main challenge we have faced in this work, is that in the classical 
linear algebra literature, b is assumed to be constant. In our settings, the 
environment is constantly changing and the computed algorithm never ter- 
minates, leading to constantly changing values of b. In this paper, we ask 
the following question: "Is it possible to devise a self-stabilizing numerical 
iterative method?" We answer affirmatively, and show that under minor 
conditions it is possible to devise a self-stabilizing algorithm that solves a 
dynamic system of linear equations, where the input to the system is con- 
stantly changing. 

We analyze the behavior of the proposed scheme under changing input se- 
quences using two different assumptions on the input behavior: a box bound, 
and a probabilistic bound. We show that in both cases, once the network 
stabilizes there is a guarantee on the quality of the computation. 

One of the most efficient distributed approaches for solving a set of linear 
equations of the type Ax = b is by using linear iterative algorithms. Unlike 
Gaussian elimination, which has a cost of 0(n 3 ), where n is the number of 
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variables, an iterative algorithm usually solves a system of linear equations 
in time of 0(n 2 r) where r is the number of iterations, which depends on the 
properties of the matrix A. In many real life problems, the convergence is 
logarithmic in n, for example in multiuser detection [5]. These algorithms are 
naturally distributed and work well in asynchronous settings. Furthermore, 
when converging, the algorithms converge to a solution from any initial state. 
An excellent survey of such methods is found in [2]. 

The main novel contribution of this paper is in analyzing the self-stabilizing 
properties of algorithms from the linear iterative methods domain. For ex- 
ample, in a practical setting, it is highly unreasonable to assume that sensor 
readings do not change over time. However, it is reasonable to assume that 
at steady state the change in sensor readings is bounded. Informally, in this 
work we show that once the input readings are bounded, the output solu- 
tion is bounded as well. This useful observation enables us to tie together 
numerical iterative methods and dynamically changing environments in a 
self-stabilizing manner. 

To the best of our knowledge, this is the first work tackling this challeng- 
ing problem. We believe that our approach can have numerous applications 
in the field of distributed self-stabilizing computation. 

Other works discuss fault tolerance aspects of distributed computation. 
For example, overcoming faults in sensors by averaging the input was investi- 
gated in [15], using a centralized algorithm. Quantifying faulty nodes' effect 
on the system's output is discussed in [11] and [8]. These papers consider 
bounded input paths and their effect on the stability of the output. In [9] 
infinite input paths are considered under the assumption that only specific 
sensors' input may change. All three papers consider discrete input values, 
as opposed to a continuous set of input values discussed in this paper. 

The paper is constructed as follows. Section 2 defines the model and 
presents the problem definition. Section 3 discusses measurement error types 
and characterization. Section 4 presents our novel algorithm SS-Iterative. 
Section 5 analyzes our algorithm and gives bounds on the convergence rate, 
assuming a box bound on the inputs. Section 6 extends our construction to 
the asynchronous case. Section 7 extends the analysis of SS-Iterative un- 
der a probabilistic model when the input is characterized using a normal dis- 
tribution. Section 8 presents experimental results of running SS-Iterative 
using sample topologies. We conclude in Section 9. 
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2. Model and Problem Definition 

Consider a distributed system of sensors measuring real-world data. Sen- 
sors are located in different areas; for example, the senors are spread through- 
out a building and they measure the temperature to adjust the heating or 
cooling. We would like the collected data to be as reliable as possible, re- 
flecting closely the changing environmental conditions. One of the obstacles 
we face when designing algorithms that collect data from a sensor network 
are measurement errors. There are two main types of inaccuracies of sensors' 
measurements: noisy environment and sensing equipment which is not cali- 
brated. In our setting sensors are allowed to communicate among themselves 
and to use data from other nodes for calibrating their environmental readings. 
Furthermore, we would like our calibration algorithm to have fault-tolerance 
properties. Specifically, we are interested in self-stabilizing algorithms [10], 
which converge to an optimal solution from any initial state. Observe that 
self-stabilization helps also in deploying the sensors. There is no need to ex- 
plicitly synchronize the sensors, once enough of them are deployed and begin 
functioning the results will converge to the expected value. 

We model the sensor calibration problem as follows. Given a directed 
communication graph G = (V,E), V is the set of nodes V = {p±, . . . ,p n }, 
E is the set of weighted edges connecting them (weights can be negative). 
Edge weights are used to model the directional dependence between nodes' 
outputs; i.e., if w PitP . = then there is no edge between and pj and their 
output is not directly dependent on each other. In addition, it is possible to 
have a non-zero self connected edge, w Pi>Pi ^ 0, which represents the weight 
of p^s own input. Note that the graph does not have to be symmetric, which 
means that it is possible that w PhP , ^ w p ^ Pi . We denote N(pi) as the set of 
node pj's neighbors, excluding p^. 

Initially, we assume a synchronous system: during a single round of com- 
munication, any pair of connected nodes may send a single message on each 
directed edge. In each round r, each node pi has a scalar input value I Pi (r), 
which represents the local reading of the sensor. 2 In addition, pi outputs its 
output value, which is denoted by Pi (r); both inputs and outputs are from 
the domain of real numbers. Denote by I(r) the input vector of the entire 
system at round r, and by 0(r) the output vector of the system at the end 



2 For simplicity of notations we use scalar variables in the paper. An extension to the 
vector case (where each sensor measures a set of measurements) is immediate. 
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of round r. In Section 6 we relax the assumption of synchronous rounds and 
provide a variant of the algorithm that works in asynchronous settings. 

The schematic operation of each node Pi at round r is composed of the 
following steps: (a) read the value of / Pl ( r )i (b) send messages; (c) receive 
messages; and (d) do some processing and output Pi {r). Then a new round 
is started, and the nodes continue so forever. 

Definition 1. A configuration C of the system at round r consists of the 
state of each node prior to performing any operation at round r; this config- 
uration is denoted by C(r). 

Definition 2. An input sequence X of length £ is a list of £ vectors such 
that each v £ X is a possible input vector of the system (i.e., v £ D, the 
domain of allowed values). An output sequence O of length £ is a list such 
that each v £ O can potentially be an output vector of the system. 

Definition 3. A step from configuration C to configuration C on input vec- 
tor v is legal if C is reached from C by the system when having v as the 
input vector, u is produced by a legal step if u is the output vector of the 
system resulting from such a legal step. 

Definition 4. A run of a system on input sequence X = {v(l), . . . , v(£)} 
starting from configuration C(r) is the sequence C(r),0(r),C(r + l),0(r + 
1),... s.t. foranyi>0: the step from C(r + i) to C(r + i + l) on input vector 
v(z + 1) is legal, and 0(r + i) is produced by that legal step. The system is 
said to produce the output sequence O = {0(r), . . . , 0(r + £ — 1)}. 

In the special case when the sensor observations (the input to the system) 
are fixed, the output decision of the sensors should converge to a solution 
that preserves the linear relations among node inputs and outputs. More 
formally, consider an input sequence X of identical input vectors; i.e., X = 
{v, v, v, . . . }. It is desired that for such an input a run from any configuration 
C on X would end up producing an output sequence O = {u(l), u(2), . . . } 
such that ||u(i) — u| | — > as i — > oo, for a u that solves the following linear 
system of equations: 

Pj&Nipi) 
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We assume that the above set of equations are uniquely solvable, denoting u 
as the solution to v. 

The following definition bounds the change in sensor observations: 

Definition 5. An input sequence X = {v(l), v(2), . . . , v(f)} is 5-bounded 

around v if for every i, 1 < % < i, it holds that ||v(z) — vj | < 5. 3 

Definition 5 states that a sequence X is 5-bounded if all the vectors in X 
are bounded within an n dimensional hypercube with an edge 25, centered 
around a point v. We note that once changes in the input are not bounded, 
then no efficient algorithm (especially in a network that is sparsely connected) 
can calculate the output fast enough. Thus it is easy to construct an example 
where the diameter of the communication graph is T>, for some system of 
equations it would take at least T> rounds for the information exchange for 
input readings at one side of the network to propagate to the other side of 
the network. 

Definition 6. Let X be an input sequence that is 5-bounded around v and 
let u be the solution to input v. A run from configuration C on input se- 
quence X e-converges to its solution if the produced output sequence O = 
{u(l), u(2), . . . u(At)} satisfies that \ \u(At) — u.) | j ^ < e(At, 5, C); where e is 
a function of At, 5 and C. 

Definition 6 requires that if - starting from configuration C - the inputs 
are in an n dimensional hypercube of radius 5 around v then the output at 
time At is bounded within some n dimensional hypercube around u with 
radius e(At, 5,C). We aim at an e(At, 5,C) that decreases as At increases, as 
long as the inputs are bounded by the same v-centered, 5-radius hypercube. 
Clearly, for 5 > 0, e(At, 5,C) > for any At. That is, there is some minimal 
radius 5' > around u and we cannot ensure a tighter bound. 

The above definition considers a single initial configuration, and a single 
input sequence X. We are interested in an algorithm that works for all initial 
configurations and all input sequences. 

Definition 7. An algorithm A e-converges for 5-bounded input sequence X 
if every run (from any configuration) on I, e-converges to its solution. An 
algorithm A e-always converges if for every 5-bounded input sequence X, 
A e-converges. 



3 H X IL = max^lxil}. 



(3 



Definition 7 formally defines the problem at hand, as an algorithm A that 
always converges has the desired self-stabilizing property: for any system 
state, once the sensors' readings changes are bounded, the change in output 
of the entire system is bounded as well. 

Our goal is to find an algorithm A that is e-always converging for a 
provably "good" e. Moreover, we aim at having A efficient also in its message 
complexity and the simplicity of the code, allowing lightweight sensors to 
actually implement it. 

3. Measurement Errors 

As mentioned in the introduction, sensor readings suffer from two types 
of inaccuracies: measurement errors and calibration errors. In this section 
we elaborate on two common types of measurement errors. This introduction 
serves as a motivation for bounding the errors using a normal distribution, a 
construction that is discussed in detail in Section 7. 

Most physical sensors tend to suffer from noise, which is defined as fluctu- 
ations in external factors to the stream of target information (signal). Most 
sensor noise can be attributed to a small number of common noise patterns, 
this noise is usually due to the nature of physical substance being measured 
and the nature of the measurement equipment. Two of the most common 
noise patterns are Shot Noise and Nyquist Noise as defined at [6]. 

Definition 8. Shot Noise is a type of noise common to the particle-like 
nature of the charge carriers. It is often thought that a DC current flow in 
any semiconductor material is constant at every instant. In fact, however, 
since current flow is made up of individual electrons and holes, it is only the 
time average flow of these charge carriers that is seen as a constant current. 
Any fluctuation in the number of charge carriers at any instance produces a 
random current change in the instance. 

As these particles (usually photons or electrons) distribute with Poisson dis- 
tribution the noise itself also follows this distribution. Poisson distribution is 
nearly Gaussian and can easily be approximated by a Normal distribution. 

Definition 9. Nyquist Noise (also Thermal noise): In any conducting 
medium whose temperature is above absolute zero (0 Kelvin), the random 
motion of charge carriers within the conductor produces random voltages and 
currents. These voltages and currents produce noise. 
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This noise is nearly white by nature (equally distributed through the spec- 
trum). It is characterized by a normal distribution of its amplitude. 

Both noise types discussed above are nearly Gaussian [12]. Therefore tol- 
erance to Gaussian noise is an important issue when examining the sensor 
calibration problem, specifically, and for systems in general. Furthermore, 
the central limit theorem states that in the average case any i.i.d sample is 
bound to be distributed as a Gaussian surrounding the mean. This makes 
Gaussian tolerance a powerful tool when dealing with dynamic input se- 
quences. 

4. Our Proposed Solution 

Getting back to the sensor calibration example, let W be the matrix that 
has Wij as entries, Eq. (1). It be written in linear algebra notation, (s.t. it 
applies to all nodes simultaneously): 

Wu = v . (2) 

If we consider a non-self-stabilizing system in which the inputs do not 
change (that is, the input is fixed, as some v), then Eq. (2) can be seen as 
Ax. = b, where A and b are given. In such a case, we are interested in 
finding the value of x, which is a vector of n unknown variables. However, 
we are interested in the case where v changes over time, and thus Eq. (2) 
does not describe the problem properly, but rather helps in understanding 
the motivation for our solution. 

Rearranging Eq. (1) we get: 

Pl (r + !) = w n>Pl I Pl {r + 1) + ^w Pi>p .O p .{r) . (3) 

Clearly, for the case of 5 = 0, a 0-bounded input sequence X, if {0 Pi {r + 
1) — Pi {r)) — > as r — > oo then Eq. (3) converges to the solution of 
Eq. (1). Thus, if the update rule of Eq. (3) is executed simultaneously by 
all nodes, and for all of the nodes (0 Pi (r + 1) — O p .(r)) — > 0, then it also 
solves Eq. (2). That is, if each node locally executes Eq. (3) then the global 
solution is reached. This observation motivates algorithm SS-Iterative in 
Figure 1. 

Remark 1. In SS-Iterative there is no notion of the "current round num- 
ber r " . That is, pi reads and writes to the variables I Pi and Pi without being 
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Algorithm SS-Iterative 



01 : Each round do: /* executed on node pi * / 

/* send current value of Pi to all neighbors * / 
02 : for each Pj G N(pi) 
03 : send Pi to pj ; 

/* update Pi according to values sent by neighbors * / 
04: set Pi := w PiiPi I Pi ; 
05 : for each value Pj received: 
06 : update Pl := Pt + w PltPj O Pj ; 

07: od. 



Figure 1: A self-stabilizing iterative algorithm. 



"aware" of r. When we discuss the algorithm "from the outside", we will 
consider I Pi (r) and Pi (r) instead of just I Pi ,0 Pi . 

Consider pi as running at some round r + 1. When pi performs Line 03, 
it sends the value of Vi . The last time Vi was updated was at Line 04 and 
Line 06 of round r. Thus, the value sent by p,i at round r + 1 is actually 
Pi (r). Therefore, the values received from pj by Pi and used to update 
Pi {r + 1) are Pj {r). However, the value read by pi in Line 04 is the value 
of I Pi {r + 1). Concluding that pi updates Pi {r + 1) exactly according to 
Eq. (3). 

Remark 2. Each node pi must know the values of w PitPj as "part of the 
code". Thus, these values cannot be subject to transient faults. 

5. Analysis of SS-Iterative 

Paper [2] shows that the update rule of Eq. (3) can be written in linear 
algebra form as 

0(r + l) = AI(r + l) + BO(r) , (4) 

where A is a diagonal matrix with w PuPi in the main diagonal, and Bij = w Pi)Pj 
for i j 

A 4 (diag(jy))- 1 , B 4 (J nxn - diag(H/)- 1 H/) , (5) 
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where I nxn is the identity matrix. Using this update rule to solve a set of 
linear equations iteratively is known as the Jacobi algorithm. 

Let X be an input sequence of length £ that is 5-bounded around vector v. 
At some round r, X = I(r), J(r + 1), . . . , I(r+£—l). Note that SS-Iterative 
saves a single scalar variable at each node, and thus the configuration of 
round r + 1 can be defined by the value of 0(r) at round r. Consider SS- 
iTERATIVE's run, starting from an arbitrary configuration at round r. We 
aim at showing that 0(r + At) is bounded by a hypercube centered at u. 
Denote by c(At) = 0(r + At) — u. If we show that ||c(At)|| is bounded 
(as At increases), then 0(r + At) is within a bounded hypercube centered 
at u. Consider c(l): 

c(l) = 0(r + l)-u 

= AI(r + 1) + BO{r) - (Av + Bu) 

= A(I(r + 1) - v) + B(0(r) - u) 

= A(I(r + 1) - v) + Bc(0) . (6) 

Since X is a 5-bounded input sequence around v, each I(r + At) can be 
denoted as v+£>(r+ At) s.t. £>(r+At) G M. n is a vector, and | \V(r + At) | < 
5. 

Claim 1. At round r + At, it holds that 

c(At) = B J AV(r + At - j) + B At c(0). 

3=0 

Proof: Proof by induction. The base of the induction was shown for 
c(l); see Eq. (6). Assume that the claim holds for At = k. Thus, c(k) = 
X^-i B j AV(r + k-j) + B k c(0). By the update rule in Eq. (4), we have that 
0(r + k + 1) = AI(r + k + 1) + BO(r + k). Combining the two equations 
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implies 



c(k + l) = 0(r + k + 1) - u 

= AI(r + jfc + 1) + BO(r + k) - (Av + flu) 

= A(I(r + fc + 1) - v) + B(0(r + jfc) - u) 

= AD(r + jfc + 1) + Bc(k) 



fe-i 



AD(r + jfc + 1) + B j+1 AV(r + jfc - j) + fi fc+1 c(0) 

3=0 
fc 

AD(r + jfc + 1) + Y B j AV(r + k+l-j) + B k+1 c(0) 



3=1 

k 

= B j AD(r + k + l-j) + B k+1 c(0) . 

Thus, if the claim holds for At = k it also holds for At = k + 1; and we have 
that the claim holds for all At > 0. □ 

Definition 10. A matrix M nXn is diagonally dominant if\Ma\ > EJ/jMjjl 
A matrix M nxn is normalized diagonally dominant (normalized, for short) 
if M is diagonally dominant, and |M„| > 1. 

Lemma 1. For a normalized diagonally dominant matrix W , it holds that 
Halloo — 1 an d loo < 1) where A,B are defined in Eq. (5) and \ \A\\ C 

\\Ax\\ 

max X7 , ^. 



I oo 



Proof: A is zero except for its main diagonal for which A^i = w 
Since \Wi_i\ > 1, we have that lA^] < 1. Thus, it holds that ||Ax|| < 



Pi-Pi 



ll x lloo- Furthermore, max x ^o TTTT^ — 1> ^- e -> Halloo — 1- Regarding B, 

I l X l loo 

Bi j = w p . iPj for i ^ j and for i = j. Since W is assumed to be normalized 
diagonally dominant, we have that \ < \W iti \, thus I^Wyl < 

1. Therefore, |-Bjj| = l w Pi,Pjl < 1 f° r a ^ I n total, for any x we 

have H-Bxlloo < H^IL, leading to H-B]]^ < 1. □ 
If W is a diagonally dominant matrix then node p^s own input effects 
Pi's output more than the sum of all of pj's neighbors outputs. That is, the 
weight of pj's input is at least the sum of weights of p^'s neighbors outputs. 
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Theorem 1. Given a normalized diagonally dominant and invertible W , 
there are constants ci,C2, where c\ > 0, and 1 > C2 > 0, such that SS- 
Iterative e-always converges with e(At,5,C) = 5ci + (c 2 ) At \ \0(r) — u^. 

Proof: By Lemma 1 it holds that II All < 1 and IIS II < 1. Consider 

•J III loo — 1 1 Moo 

a 5-bounded input sequence X around v, and SS-Iterative's run starting 
from an arbitrary state 0(r). We are interested in the behavior of | |c( At) \ \ : 



\c(At)\\ 



< 



< 



< 



B At c(0)\ 



At-l 

B j AV{r + At- j) + B At c(0) 

3=0 
At-l 

B j AV(r + At — j 

At-l 

E WBWiWAVir + At - j) 

j=0 



B 



I At 



|c(0)|| c 



5\\A\ 



S\\A\ 



At-l 

i y 

loo / j 
3=0 

1 - | 

loo 



\B\ 



\B 



i At 



LB 



\B 



\B 



i At 



i At 



|c(0)|| c 



(7) 



For an input sequence X that is 5-bounded around v, denote by u the solution 
to the original system of equations Wu = v. By Eq. (7), 

|0(r + At)-u|| <5|L4|I 



1 

J 



Since H-B^ < 1, we have that !_ B °° < j 

I At 



LB 



i 



|iC*||c(0)|l 



it holds that \\A 
c(0) = 0(r 



l- 



oo 1-||B,_ 

u we are done. 



and by setting ci 

. . II oo 

< c\. By setting c 2 = ||£> 



Ml 



Halloo 

and recalling that 

□ 



Theorem 1 states sufficient conditions s.t. SS-Iterative e-always con- 
verges. Moreover, the algorithm SS- Iterative is lightweight, as it requires 
nodes to send only a single value to every neighbor on each round. 



6. Extension to the Asynchronous Model 

Our second novel contribution is in extending our model to support asyn- 
chronous communications. In a large sensor network, it is unreasonable to as- 
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sume that the sensors operate in synchronous rounds. Furthermore, as known 
from the linear iterative algorithms literature, algorithms sometimes converge 
in less asynchronous rounds (when compared to synchronous rounds). 

When considering the asynchronous model, it is more convenient to dis- 
cuss shared-memory as means of communication. 4 Thus, assume that for 
each directed edge between Pi,Pj there is a read- write register R PuPj that is 
written by pi and read by pj . 

An asynchronous run is an infinite sequence of configurations Cq — > C\ — > 
. . . such that some process p performs an atomic step between configuration 
Ci and Ci+i. An atomic step consists of reading or writing from a single 
register. Notice that in the current model a configuration consists of all the 
registers and of the local variables at the different nodes. 

In this section we again prove that starting from an arbitrary configura- 
tion, when the changes to the inputs are bounded, the outputs are bounded 
as well. We consider each configuration C r to be assigned a vector input I(r) 
such that if node pi reads the input when performing an atomic step on C r 
it reads the value of I Pi {r). Equivalently, the output vector of configuration 
C r is O(r). 

Figure 2 outlines Async-SS-Iterative, which is a direct translation of 
SS-Iterative to the shared-memory model. 

Async-SS-Iterative consists of two phases: in the first, the previous 
value of Pi is written to all its neighbors' registers. In the second phase pi 
calculates its new value of Pi by reading the registers of all its neighbors. 

We consider only "fair" runs, in which each node performs an atomic step 
infinitely many times. Thus, each node performs both phases infinitely many 
times. A round is defined to be the shortest prefix of a run such that each 
node has performed at least once lines 02-07 in the algorithm. We number 
each atomic step and each round. Note that a round consists of many atomic 
steps. 

We model a fair run as follows. Each node Pi performs infinitely many 
atomic steps, and participates in infinitely many rounds. Notice that the 
registers pi reads in round k + 1 have all been last written to, no earlier than 
during round k. Since a round consists of each node performing all the steps 
in the algorithm, each node Pi manages to read all of its neighboring registers 



4 In [10] it is shown how to convert an algorithm based on shared-memory to a message- 
passing algorithm with links of bounded capacity. 
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Algorithm Async-SS-Iterative 



01 : Forever do: 



/* executed on node pi */ 



02: 
03: 



/* write current value of Pi to all neighbors 
for each Pj e N(pi): 
write (J Pi to R Pi:Pj ; 
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04 
05 
06 
07 



/* update Pi according to values of neighbors */ 
for each pj e N(pi): 



read R Pj . Pi into temp; 



update Pi := Pi + w puPj temp; 



08: od. 



Figure 2: A self-stabilizing iterative algorithm for asynchronous networks. 



and to write to all of its neighboring registers every round. Thus, there is 
some atomic step r (during round k + 1) such that: 



where r', r'j (for all pj ^ pi) are smaller than r and are from at least round k. 

Let u be such that u = Av + Bu, and let the inputs be from a 5-bounded 
input sequence around v. Denote c(r) = 0(r) — u and z = maxj |c p .(0)|. 

Theorem 2. Given a normalized diagonally dominant and invertible W , 
and while considering only fair runs, there are constants c±, C2, where c\ > 0, 
and 1 > C2 > 0, such that Async-SS-Iterative e-always converges with 
e(At,5,C) = 5ci + (c 2 ) At z; where At counts the asynchronous rounds of a 
fair run. 

Proof: Notice that if pi did not perform the rth atomic step then Pi (r) = 
Pi (r — 1) and therefore c Pi (r) = c Pi (r — 1). Consider the value of c Pi (r) 
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when pi did perform the rth atomic step (during round k + 1). 



c vA r ) = Pi (r) 



w Pi ,p 



W Pi,Pj U Pi 



' P A r j)- u Pi 



!>, C !',( r j) 1 



where r' and the different r'- are smaller than r and are all from round k or 
round k + 1. 

By using Lemma 1 we get: 

|c Pl (r)| < |w PiiPi (/ Pi (r') - v Pi )| + max |c p .(rp| V |w Pi)P .| 

Pj ... 

< \ W p i ,Pi{ I Pi{ r ') ~ Y Pi)\ + Halloo \ C Pma*( r max)\ 
— ^ "f" I l-^l loo l C Pmai ( r maz) I ) 

for some p max . and r max < r that is from round k or k + 1. 

Therefore, for any p 4 during round fc + 1 there is a list of length £ > k 
of nodes Pi,P2, ■ ■ ■ ,Pe an d a sequence of length £ of atomic steps T\ > r 2 > 
• ■ ■ > re — 0, such that 

|c Pi (r)| < S+WBlUc^n)] 

< *+l|S|| 00 (*+||S|| 00 |c B (r 2 )|) 
= S(l + \\B\U + \\B\\l\c P2 (r 2 )\ 



< SJ2\\B\\l + \\B\t\c n (n)\ 

2=0 

a - \\b 



2 = 

|£CM0)|. 



I DC 



ll^ll 



oo 



Denote by c\ = impm , and c 2 = H-B^. We have that for node pi 

II M oc 

performing the rth atomic step during round k it holds that |c Pi (r)| < 8c\ + 
c^z < 5ci + c 2 z. □ 
In fair runs, there are infinitely many rounds k, thus, as / and r go to 
infinity, we have that ||0(?")lloo ^ s bounded by a hypercube of length 8c\ 
around u. 



15 



7. Probabilistic Bound 

In this section we analyze the behavior of the SS-Iterative algorithm, 
when the input readings are normally distributed due to noisy readings. 

Let X be an input sequence of length I that is normally distributed around 
a mean vector v. That is, X = I(r), J(r + 1), . . . , J(r + I — 1), where each 
is sampled i.i.d from a Gaussian distribution with mean v and covariance 
S v . Consider SS-Iterative's run, starting from an arbitrary configuration 
at round r. We aim at showing that 0(r) is equivalent to a sequence O 
sampled i.i.d from a Gaussian distribution O ~ A/"(u, S u ) where u is the 
solution to v. 

Since we know the distribution X, we can calculate 0{r) using Eq. (4) 
after I{r) is sampled. O(r), our sequence's starting point, can take arbi- 
trary values as we will show that SS-Iterative converges from any starting 
position. 

We can now deduce that after k rounds of sampling the input distribution 
the output 0{k + 1) is distributed as: 

0(k + 1) = AI(k + 1) + BO{k) ~ Af(Av + BO{k), AE V A T ) . (8) 

This formulation is, however, dependant upon O(k), or upon the specific 
samples of /(•). The output of each round is normally distributed with mean 
Av + BO{k) and with covariance matrix AT, V A T . 

To find O we examine the recursive formula for 0(k+l). For simplifying 
notations, we mark the starting round r as round 0. 

0{k+l) = AI{k + l) + BO{k) 

= AI{k + l)+BAI{k) + B 2 0{k-l) 

= AI{k + 1) + BAI{k) + ... + B k {AI(l) + BO(0)) . 
Now we take k — > oo to find the distribution O: 

O oc lim (I nxn + B + B 2 + ... + B k )AX + lim B k+1 O(0) . (9) 

k— »oo k— *oo 

By Lemma 1 it holds that < 1, lirn^oo B k = 0. Using the Maclaurin 

series we get: 

lim I nxn + B + B 2 + . . . + B k = (I nxn - B)-\ 

k— >oo 
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Denoting B = (I nxn ~ B) 1 we write the distribution of O: 

O oc BAT ~ Af(BAv, AB T J: V BA) . 

Recall the definitions of A and B (5) and find that: 

BA = (I nxn - By 1 A = (I nxn + AW- I nxn )- l A = (AW)" 1 A = W~ l . 

Note that the output sequence can now be characterized more intuitively as: 

O oc W~ X X ~ Af(u, W r - 1 E v W -lT ) . 

Using known bounds such as multivariate extensions of Chebyshev's bound 
one can now introduce bounds similar in nature to the box bounds discussed 
earlier. A better bound can be acquired using the cumulative distribution 
function for normal distributions; however these can only be estimated nu- 
merically. The use of such bounds is: 

p = Pr(||/(fc) - v|| < 5) , 3e\ Pr(||0(fc) - u|| <e)=p t 

where the connection between 5 and e depends on the specific bound used. 

So far we have shown that SS-Iterative is tolerant towards Gaussian 
noise, and the output sequence converges to a sequence drawn i.i.d from a 
normal distribution around the expected solution. 

7.1. The Asynchronous Setting 

Our proof so far discussed the synchronous case. Next, we extend the 
analysis of the probabilistic bound for asynchronous case, under mild as- 
sumptions. We assume the "Totally Asynchronous" model as discussed in 
[3]. We denote T l as the set of times in which the i'th node updates its 
neighbors with new output values. We further denote T % {tk) as the vector 
containing the last update times of all received messages at node i at time 
tk- 

Definition 11. Total Asynchronism: the sets T l are infinite and if {tk} 
is a sequence of elements ofT 1 that tends to infinity, then linifc-joo rj(tfc) = oo 
for every j . 
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This assumption implies that every node is updated infinitely often while 
allowing maximal flexibility for delay and loss of data. 

[3] gives a converges proof for the asynchronous case, to an iterative al- 
gorithm with constant input values. The challenge we face is that the input 
is not constant, but distributed normally. Thus we need to prove that the 
solution does not infinitely diverge. 

For proving that the algorithm output does not diverge, we require some 
further assumptions. We assume that the update times of each of the nodes 
are independent of the distribution X. The second assumption we make is 
that the input noise distribution is normally distributed, where the noisy 
reading in different sensors are not correlated. In other words, the inverse 
covariance matrix of the input is a diagonal matrix. 

To prove that under these assumptions O is normally distributed we 
define the set T: 

T = {t\3i:teT}. 

Denote = 0, ... as the ordered sequence of times in T. For all G T we 
need to find the distribution of the input at time t^. 

Theorem 3. The infinite input sequence I(tu)) is distributed normally with 
X. 

Proof: As X is normally distributed, any marginal distribution computed 
by a subset of the nodes is normally distributed as well. Now we can induce 
the requested theorem: 

Induction base: By the Total Asynchronism assumption there exists a 
time to > at which all nodes have updated their value, and thus the input 
of each node is now drawn from its marginal distribution. This means that 
at time to the distribution of the entire input vector is (by cartesian product) 
drawn from X. 

Induction hypothesis: At time t = t^ we assume the input distribution 
is distributed normally from X. 

Induction step: At time t = we separate the nodes into two groups: 
The first group contains nodes that remained constant from the previous 
round, that is, their a-priori distribution is their marginal distribution drawn 
from X by our induction hypothesis. 

The second group contains nodes that are updated with a new input value 
at this round; this value is distributed again from their marginal distribution. 
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Notice that all of these marginal distributions are Gaussian, thus, their 
cartesian product is Gaussian as well. In addition, any joint covariance be- 
tween nodes is due to our assumption. Therefore, we have a cartesian 
product of all of the nodes distributed as marginal distribution of X and thus 
the entire input I(t^ + i)) follows the distribution X. 

Finally we induce that at all tn) G T: /(£(;)) are drawn from X. Note that 
unlike the synchronous model - these samples are no longer taken i.i.d from 
the distribution X. □ 

Now we can use Theorem 3 and the Total Asynchronism assumption to 
draw a similar proof to the synchronous case. Recalling Eq. (9): 

O oc lim (I nxn + B + B 2 + . . . + B k )AI + lim B k+1 O(0) . 

k— >oo k— *oo 

Theorem 3 shows that the input distribution, at all rounds, is drawn from X. 
This implies that we may apply the Maclaurin series again: 

lim I nxn + B + B 2 + . . . + B k = (J nxn - B)-\ 

k~*oo 

We now turn to [3, Section 6.3.2] where the following claim is proven: 

Lemma 2. For a normalized diagonally dominant matrix ^A, it holds that 
x = (I nX n — , yA)x + 76 asynchronously converges. 

We now make use of Lemma 2 under the following notations: 6 = 0, which 
means that our solution converges to 0. Recalling Eq. (5): B — I nXn — AW, 
where W is diagonally dominant and A is a diagonal matrix, results in AW 
being diagonally dominant. Combining these results together, we conclude 
that B is a contraction, thus the linear dynamical system x = Bx converges 
to zero for any initial x: 

lim B k+l = . 

k— *oo 

Regarding convergence speed, the speed at which the initial output factor 
converges to and the output distribution converges to a normal distribu- 
tion, is dependant upon the rate of the least updated node. These settings, 
though somewhat more restricting than the ones used for synchronous con- 
vergence, are all "minimal" in the sense that removing any of them creates a 
system that may diverge. Further discussion about the validity of the above 
assumptions is given in Section 9. 
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8. Experimental Results 

8.1. Box Bound 

For illustrating the behavior of our proposed algorithm using the box 
bound assumption, we have simulated SS-Iterative using two sample topolo- 
gies of one hundred nodes. Figure 3 depicts a circular topology where each 
node is connected to its left and right neighbors. Figure 4 shows a random 
unit disc graph, where nodes are randomly spread on a plane, and each node 
is connected to the nodes that are within a distance of 1. The X-axis shows 
the number of iterations, and the Y-axis shows the value of 5. Area colors in 
the heatmap depict the average of the following procedure: randomly select 
a vector v and a 5-bounded sequence around v, run the simulation for the 
randomly selected values and return the distance between the last output 
vector and u (calculated as u = W -1 v). The heatmap uses a log log scale. 
Both graphs clearly show that as 5 decreases and the number of iterations 
increases, the output of SS-Iterative converges to be bounded by a small 
hypercube around u. 

Note that the unit disc weighted topology matrix is characterized by 
Halloo = 0.02, | \B\ \ 00 = 0.97 while the circle graph is characterized by 
1 1^1 loo = 0-33, loo = 0.66. As expected, using unit disc topology requires a 
larger number of iterations for convergence (depends on ). In addition, 

in the unit disc topology the value of 5 has a smaller effect on the conver- 
gence, due to the value of H^H^, which affects the minimal radius around 
the output. Since H-A]^ is smaller in the unit disc topology, increasing S 
does not significantly affect the convergence. 
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8.2. Probabilistic Bound 

We have simulated a run of SS-Iterative, when the input noise is i.i.d. 
normally distributed. Figure 5 plots the behavior of the output distribution 
at specific rounds. The input distribution of two nodes, printed in blue, is 
a normal distribution where the two nodes are uncorrelated. The output 
(in red, growing darker in later rounds) are distributed around some round 
mean, with the means themselves distributed around u. It is important to 
note there is no greedy convergence to u. This is because each round contains 
a different noise sample and thus at a specific round our solution can stray 
away from u. 

Figure 6 plots the behavior of the mean value of the input and output to 
the large system limit. The input is to the right (in blue axes) the output 
is to the left (gray circles). After a large number of rounds, the output is 
normally distributed around u, where the plotted u (in green) is the result 
of an EM estimation of the sample mean. 

Figure 7 plots a similar experiment from a different angle. This figure 
shows the input reading of the nodes (to the right) which is spherically dis- 
tributed in a circle because the input noise is uncorrelated. The output 
of the SS-Iterative algorithm (to the left) is correlated among nodes, as 
seen by the elliptic shape of the distribution. The correlation of the output 
distribution is supported by the theoretical results. 




Figure 5: Output of several rounds. Figure 6: Convergence behavior. 

9. Discussion 

We have shown that the algorithm SS-Iterative is a modification of 
the Jacobi iterative method to solve a set of equations Ax = b, where A 
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Figure 7: Self correlation of input and output. 

is given and b is dynamically changing but bounded or statistical. More- 
over, Theorem 1 is a generalization of previous analysis of Jacobi's conver- 
gence. Our motivation for SS-Iterative originates from the sensor calibra- 
tion problem, where sensors need to calibrate their noisy readings. Unlike 
previous approaches to this problem, we assume a dynamic system with an 
infinite execution of the algorithm. In this setting the readings of the sensors 
continuously change. Under the assumption that the readings' changes are 
bounded, we have shown that the calibrated output is bounded as well. 

Further application for SS-Iterative can be found in any setting where 
it is desired to solve Ak = b in a converging and self-stabilizing manner, 
while A is given, and b may change slightly from one round to the next. 
Notice that the analysis given in Section 5 holds in such a system. 

As noted in Remark 2 the matrix A is "part of the code" . An optional 
alternative to the current solution is to compute A~ l (the inverted matrix 
of A) beforehand and include it "as part of the code". Thus, each node 
could locally solve x = A _1 h, and it can be shown that x will be bounded 
(as long as b is bounded). The main problem with such a solution is the 
connectivity requirements it incurs. In our solution, scalar values are sent 
in the network only between direct neighbors. The matrix W represents a 
weighted adjacency graph. Once inverted, the matrix W~ l might not be 
sparse. A non-zero entry w^ 1 G W" 1 means that node Pi needs to communi- 
cate with node pj. This extra communication might cause the algorithm to 
lose its self-stabilizing properties, as non-neighboring nodes would require a 
self- stabilizing overlay network for their communication. 

The assumption of a predefined A is suitable for static networks in which 
the communication graph is predetermined. For dynamic networks, it would 
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be interesting to adjust SS- Iterative to discover the connectivity of the 
network, inferring the optimal weights dynamically. We assume that after 
the weights are calculated, the topology of the sensor network remains stable, 
thus the convergence analysis of Section 5 should hold. 

9.1. Necessity of Convergence Conditions for the Probabilistic Bound 

Assuming that the sensor input readings are normally distributed due 
to measurement noise, the resulting output is also normally distributed. As 
discussed in Section 7, this is not a sufficient condition for the convergence 
of SS-Iterative in a finite number of rounds. 

We give the following example to demonstrate the claim above. Assume 
the input reading at round 1 is larger than v, and at each of the following 
rounds the input reading of each node is larger than at the previous round. 
As the input is normally distributed, for every n > there is a small but 
positive chance for such a sequence of length n. Input sequence behaving 
as described would lead SS-Iterative to a growing output sequence for n 
rounds, meaning for the duration of these rounds SS-Iterative diverges. 
Thus it is obvious the mean can be infinitely far from u as long as n is finite. 
However, when n is large, the probably of such an anomaly is diminishing to 
zero. 

The convergence results for the asynchronous case had a minimal set 
of assumptions. One of them is that the input noise distribution must be 
uncorrelated among nodes. That is because if the input noise distribution 
is correlated, all nodes must sample the input concurrently, an assumption 
which is not valid under the asynchronous model. 

9.2. Relation to Perturbation Theory 

A large amount of research focused on the problem of solving Ax = b 
when A and b are not exactly known. That is, let A = A+5A and b = b+<5b, 
and consider the equation Ax = b; what can be said about x in relation to 
x? 

Our setting is "easier" in one sense, and "harder" in a different sense. In 
our setting A is known, i.e., 5 A = 0. However, b is not well defined. That 
is, the input vector - which is described by b - changes over time. When 
solving Ax = b it is assumed that there is some b that is constant but it 
was measured with an error. In our case, b is not constant as it changes over 
time, while it represents the measurements correctly. 
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As a future research, it would be interesting to consider the implications 
of adding inaccuracy to the measurements. The vast body of knowledge 
regarding perturbation theory would definitely aid in this extension to our 
model. 

9.3. Relation to Convex Optimization 

Many practical optimization problems are given in the quadratic form 
fix) = l/2xAx — b T x, where the task is to compute min x /(x) distributively 
over a communication network. A survey showing several applications can be 
found in [4]. Example applications are monitoring, distributed computation 
of trust and ranking of nodes and data items. 

A standard way for solving min x /(x) is by computing the derivative and 
comparing it to zero to get the global optimum. When the matrix A is 
symmetric, /'( x ) = ^ x — b = 0, and we get a linear system of equations 
Ax = b. In other words, the convex optimization problem is reduced into a 
solution of a linear system of equations. 

Interior point methods [7, Ch. 11] solve linear programming problems 
by applying Newton method iteratively. Each computation of the Newton 
step involves a solution of a linear systems of equations. An area of future 
work is to examine the applications of our self-stabilizing algorithm to these 
methods. The difficulties arise from the fact that the matrix A needs to be 
recomputed between iterations, so nodes need to be synchronized and aware 
of the current iteration taking place. 
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